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ABSTRACT 

The evolution of marginally bound supercluster-like objects in an accelerating ACDM 
Universe is followed, by means of cosmological simulations, from the present time to 
an expansion factor a = 100. The objects are iden tified on the basis of the binding 
density criterion introduced bv lDiinner et al.l (|2006D . Superclusters are identified with 
the ones whose mass M > lQ^^h~^MQ, the most massive one with M ^ 8x lG^^ft,~^Mo, 
comparable to the Shapley supercluster. The spatial distribution of the superclusters 
remains essentially the same after the present epoch, reflecting the halting growth of 
the Cosmic Web as A gets to dominate the expansion of the Universe. The same trend 
can be seen in the stagnation of the development of the mass function of virialized 
haloes and bound objects. The situation is considerably different when looking at the 
internal evolution, quantified in terms of their shape, compactness and density profile, 
and substructure in terms of their multiplicity function. We find a continuing evolution 
from a wide range of triaxial shapes at a = 1 to almost perfect spherical shapes at 
a — 100. We also find a systematic trend towards a higher concentration. Meanwhile, 
we see their substructure gradually disappearing, as the surrounding subclumps fall 
in and merge to form one coherent, virialized system. 

Key words: Cosmology: theory - large-scale structure of Universe - galaxies: clusters: 
general 



1 INTRODUCTION 

The evidence for an accelerated expansion of the Universe 
has established the dominant presence of a 'dark energy' 
component. In the present cosmological paradigm, the Uni- 
verse entered into an accelerating phase at z ~ 0.7. Obser- 
vational evidence points towards a dark energy component 
which behaves like Einstein's cosmological constant. As long 
as the matter density in the Universe dominated over that 
of dark energy, the gravitational growth of matter concen- 
trations resulted in the emergence of ever larger structures. 
Once dark energy came to dominate the dynamics of the 
Universe and the Universe got into a ccelerated expansion, 
structure formation came to a halt (|Peeblesl [l98ol : iHeathI 
1 19771 ). With the present-day Universe having reached this 
stage, the largest identifiable objects that will ever populate 
our Universe may be the ones that we observe in the pro- 
cess of formation at the present cosmological time. While no 
larger objects will emerge, these sufficiently overdense and 



bound patches will not be much affected by the global cosmic 
acceleration. They will remain bound and evolve as if they 
are island universes: they turn into isolated ev o lving regions 
JChiueh fc Hel I2OO2I: iNagamine fc Loebl l2003l : iBusha et all 
bOOAlDiinner at aLfeooef T 

While clusters of galaxies are the most massive and 
most recently collapsed and virialized structures, the 
present-day superclusters are arguably the largest bound 
but not yet fully evolved objects in our Universe. In our 
accelerating Universe, we may assume they are the objects 
that ultimately will turn into island universes. A large range 
of observational studies, mostly based on optically or X- 
ray selected samples, show that clusters are strongly clus- 
tered and group together in large supercluster complexes 
(see e.g.fOort 1983; B ahcaU 1988; Ein asto et al..,1994. . ,20oTI : 
lOuintana et al.lboool '). These superclusters, the largest struc- 
tures identifiable in the present Universe, are enormous 
structures comprising a few up to dozens of rich clusters 
of galaxies, a large number of more modestly sized clumps. 



2 Araya-Melo et al. 



and thousands of galaxies spread between these density con- 
centrations. 

In this study, we aim at contrasting the large-scale evo- 
lution of structure in an accelerated Universe with that of 
the internal evolution of bound objects. In order to infer 
what will be the largest bound regions in our Universe, 
the island universes, we study the mass function of bound 
objects. The abundance or mass function of superclusters 
serves as a good indicator of the growth of structure of a 
cosmological model. While large-scale structure formation 
comes to a halt, this will manifest itself in the asymptotic 
behaviour of the supercluster mass function. Meanwhile, the 
internal evolution of the superclusters continues as they con- 
tract and collapse into the largest virialized entities the Uni- 
verse will ever contain. 

We address three aspects of the continuing internal 
evolution of bound regions: their shape, density profile 
and internal substructure in terms of their cluster mul- 
tiplicity. The shape of supercluster regions is one of the 
most sensitive probes of their evolutionary stage. We know 
that superclusters in the present-day Universe are mostly 
flattened or elongated structures, usually identified with 
the most promin ent filaments and s heets in the galaxy 
distribution (e.g. [Plionis et al.' '1992'; ^Sathvaprakash e t al.l 
19981: iBasilakos et al. 2001; Shcth, Salmi & Shandarin 200^; 



Einasto et al.l 2007c ). The Pisces- Perseus supercluster chain 



is a particularly well-k nown example of a strongly elon- 
gated filament (see e.g. iHavnes et al.l[l98^ ). The distribu- 
tion of shapes of bound structures is a combination of at 
least two factors. One is the shape of the proto-supercluster 
in the initial density field. The second factor is the evo- 
lutionary state of the bound structure. We know that the 
gravitational collapse of cosmic overdensities — whose pro- 
genitors in the primordial density perturbation field will 
never be spherical (Peacock & Heavens 1985; Barde en et al.l 
1 19861 ) — proceeds in a distinctly anisotropic fashion via fiat- 
tened and elongated configurations towar ds a final more 
compact t ri axial virialized stat e (see eg. ZeldovichI 



. 1973 : 

Icke '1973^; 'White & Silk' '1979"; 'Eisenstein & Loeb' '1993; 
Bonc L jg Myers 1996j Satliyapra . kash et al. 1996; Dosiacqueg 
Ivan de WevgaerTfc Bondl [20081 1. 
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The collapse of the superclusters will also result in a 
continuous sharpening of the internal mass distribution, re- 
fiected in the steepening of their density profile. While they 
are in the process of collapse, internal substructure of con- 
stituent clusters remains recognizable. While the subclumps 
merge into an ever more massive central concentration the 
supercluster substructure gradually fades, resulting in an in- 
creasingly uniform mass distribution. The evolving and de- 
creasing level of substructure will be followed in terms of the 
evolving supercluster multiplicity function, i.e. the number 
of cluster-sized clumps within the supercluster region. 

Representing moderate density enhancements on scale 
of tens of Mpc, in the present Universe superclusters are 
still expanding with the Hubble flow, although at a slightly 
decelerated rate, or have just started contracting. Because 
these structures have not yet fully formed, virialized and 
clearly separated from each other, it is difficult to iden- 
tify them unambiguously. In most studies, superclusters 
have been defined by more or less arbitrary criteria, mostly 
on the basis of a grouping and /or percolation algorithrn 
(see e.g IOortlll983l : iBahcalllTigSSI : lEinasto et aLlll994l . I2001I : 



iQuintana erallbOOOl ). This introduces the need for a user- 
specified percolation radius. iDiinner et al.l (|2006l . hereafter 
Paper I) attempted to define a more physically based crite- 
rion, identifying superclusters with the biggest gravitation- 
ally bound structures that will be able to form in our Uni- 
verse. On the basis of this, they worked out a lower density 
limit for gravitationally bound structures. This limit is based 
on the density contrast that a spherical shell needs to enclose 
to remain bound to a spherically symmetric overdensity. 

We use this spherical density criterion to identify bound 
structures in a large cosmological box. In this we fol- 
low the work of _^iuch & He (200^ and |d iinner et al.l 
(|2006l l. IChiueh fc Hd (20021 numerically solved the spher- 
ical collapse model equations for self-consistent growing 
mode perturbations in order to obtain a theoretical criterion 
for the mean density enclosed in the outer gravitationally 
bo und shell. The resulting density criterion was evaluated 
by iDiinner et al.l (|2006l ) on the basis of numerical simula- 
tions. They generalized it by deriving the analytical solu- 
ti on which also f orms the basis of the current study, and 
m IDiinner et al] l|2007l ) extended the criterion to limits for 
bound structures in redshift space. 

Various authors ha ve addressed the future evolu- 
tion o f cosmic structure (IChiueh fc Hel 20021; Busha et ajj 
20031: iNagamine fc Loebl l2003l: IDiinner et all [2OO6, 200l 



Hoffman et all 120071 : iBusha et all 120071 ) . The internal evo 



lution of the density and velocity structures of bound ob- 
j ects w as followed bv lBusha et al.l l|2003l ). with lBusha et al.l 
iiooj) focusing on the effects of small-scale structure on the 
form ation of dark matt e r halo s in two different cosmolo - 
gies. iNagamine fc Loebl (|2003D and iHoffman etall (|2007D 
specifically focused on th e evolution of the Local Universe. 
INagamine fc Loebl l|2003t ) found that the Local Group will 
get detached from the rest of the Universe, and that its 
physi cal distance to other systems will increase exponen- 
tially. iHoffman et al.l (|2007t ) investigated the dependence on 
dark matter and dark energy by contrasting ACDM and 
OCDM models. They concluded that the evolution of struc- 
ture in comoving coordinates at long term is determined 
mainly by the matter density ra ther than by the dark en - 
ergy. A key po i nt of atten tion in Nagamine fc Loebl l|2003l ). 
IHoffman etal] (120071 1 and lBusha et all (|2007l l" was the mass 
function of objects in their simulations, on which they all 
agree that it hardly changes after the current cosmic epoch. 

This paper is the first in a series addressing the future 
evolution of structure in FRW Universes. In an accompany- 
ing publication we will specifically look at the infiuence of 
dark matter and the cosmological constant on the emerging 
(8uper)clust e r p opulation, based on the work described in 
lArava-Meld (|2008l ) . The present study concentrates on the 
details of this evolution in a standard flat ACDM Universe. 

This paper is organized as follows. In section (2] we 
present a review of the spherical collapse model, including 
a derivation of the critical overdensity for a structure to 
remain bound. Section [3] describes the simulation and the 
group finder algorithm that we employ when determining 
the mass functions. This is followed in sect. |3] by a a qual- 
itative description of the evolution, including a case study 
of the evolution of some typical bound mass clumps from 
the present epoch to a = 100. The lack of evolution in their 
spatial distribution in the same time interval is studied in 
sect.O Section |6] presents the mass functions of the bound 
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structures at a = 1 and a = 100 and a comparison with the 
ones obtained by the Press-Schechter formahsm and its vari- 
ants. The evolution of the shapes of the structures is studied 
in section [T] In section [8] we look into the mass distribution 
and density profiles. Section |9] presents the stark changes in 
the supercluster multiplicity function. In sect. [TD] its results 
are combined with those on the supercluster mass functions 
obtained in section [6] to relate our findings to the presence 
and abundance of monster supercluster complexes like the 
Shapley and Horologium-Reticulum supercluster. Finally, in 
section [TTl we discuss our findings and draw conclusions on 
various issues addressed by our study. 



2 SUPERCLUSTER COLLAPSE MODEL 

The present study is based on a physic al criterion for th e 
definition of superclusters proposed bv iDiinner et al.l (|2006l ). 
namely that they are the largest bound (though not yet viri- 
alized) structures in the Universe. The practical implemen- 
tation of this definition suggested in that paper is through an 
(approximate) density threshold for regions that are gravita- 
tionally bound. Given the anisotropic nature of t he collapse 
and t h e tidal influ e nce of the cosmic surro u ndings (|Zeldovich| 
1970| : llckd 1 1971 'Bond fc Mver^ 1 19961: ISheth fc TormenI 
19991 : iDesiacaueslB oOS: van de Weygaert fc Bond 2008) this 
may only yield a rou gh approximation. The work by 
ISheth fc TormenI l|2002i l did show that a density threshold 
does depend on shape and environment. One could imag- 
ine a variety of alternative physical definitions for super- 
clusters, also relating to the assumption that they are the 
largest bound structures in the Universe which just have 
commenced to condense out of the cosmic background. One 
particular criterion would be to invoke the corresponding ve- 
locity field and identify them with bound regions that have 
turned around and started to contract, on the way towards 
complete collapse a nd virialization. 

Here we follow iDiinner et all l|2006l ) and assume that a 
global density threshold criterion, in combination with a few 
extra assumptions, assures a reasonably accurate identifica- 
tion process. This has indeed been demonstrated in the same 
study in a comparison of the criterion with the outcome of 
nu merical simulat io ns. 

iDiinner et all l|2006h derive an analytical density 
threshold criterion for bound regions in a Universe with 
dark energy. We summarize this criterion and derivation in 
the subsections below. In addition to the density threshold, 
we assure that a given bound region has started to mate- 
rialize as a recognizable entity by including the additional 
requirement of the bound regions to have a virialized core. 
In a final step we group the overlapping identified spheri- 
cal bound supercluster objects in order to outline a region 
that in the observational reality would be recognized as a 
supercluster. 

The spherical density criterion described below forms a 
key ingredient of our supercluster definition and for the iden- 
tification procedure for singling out bound spherical regions 
associated with clusters (sect. [37l.5[) . 



2.1 Spherical Collapse Model 

The spherical coUapse model l|Gunn fc GottI Il972l : 
iLilie fc Lahavlll99ll : Eahav et al.lll99lh describes the evolu- 
tion of a spherically symmetric mas s density p erturbation 
in an expanding Universe (also see iMota fcva n do Bruc^ 
|2004 ). Its great virtue is the ability to completely follow the 
nonlinear evolution of a collapsing shell, as the dynamics is 
fully and solely determined by the (constant) mass interior 
to the shell. Even though the gravitational collapse of 
generic cosmological structures tends to be highly inhomo- 
geneous and anisotropic, the spherical model has proven 
to provide a surprisingly accurate description of the more 
complex reality. Even in situations where it is not able to 
provide accurate quantitative predictions it may be used as 
a good reference for interpretation of results. 

We consider a flat FRW Universe with a cosmologi- 
cal constant A. The mass density parameter at the current 
epoch is Qm,o, while nA,o is the present value of the cos- 
mological density parameter. The Hubble parameter at the 
current epoch is Ho- 

A mass shell with a (physical) radius r{t) at time t 
encloses a fixed mass M. The starting point of our derivation 
is the energ y E pe r unit mass of the shell, which satisfies the 
equation (lPeebQ[l983 ): 



2\dt 
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: constant . 



By introducing the dimensionless variables 
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the energy equation can be recast into the simpler dimen- 
sionless form 



E^(9^\ E. 
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(5) 



2.2 Critical shell and turnaround radius 

To delineate a bound regions around a spherical mass con- 
centration we have to identify the critical shell. It is the shell 
which separates the regions that will expand forever and the 
ones that will at some moment in time turn around and fall 
in onto the core of the region. 

In a Universe with a cosmological constant, the critical 
shell is the one that delimits the region of gravitational at- 
traction and the region of repulsion. This translates into the 
radius for which the (dimensionless) potential energy V, 



r 2 



(6) 



is maximized. The maximum of this potential occurs at f* — 
1. The critical shell (indicated by the subscript 'cs') is the 
shell with the maximum possible energy to remain attached 
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to the spherical mass concentration, 

= y(f* = l) = (7) 

The evolution of the critical shell's radius fcs(t) can be in- 
ferred by integrating the energy equation ([5]) from t — 
(f = 0) until epoch i, for the energy E = E* — of the 
critical shell. In a flat Universe with a cosmological constant, 
the vacuum energy density parameter SIa is a monot onically 
increeising function of the age t of the Universe l|Peeblesl 



, 2 / 3t 
tanh — 



(8) 



The evolution of the critical shell fcs(t) may therefore be 
expressed in terms of i}A{t), 



Xir) 



_xif) + i 

where the variable xi^^) is given by 



(9) 



X(r) 



1 + 2r + ^3f (f + 2) 
1 + 2f - ^3f (f + 2) 



(10) 

Evaluation of this expression shows that the critical shell 
radius fcs is an increasing function of Qa- The shell reaches 
its maximum - turnaround - radius r-max as Qa ^ 1 and 
t — » oo. The maximum radius is reached when 1/(fmax) = 
E* — — |, i.e. when 



(11) 



so that the normalized radius can be interpreted as 

For the cosmology at hand {Q,m.o ~ 0.3, SIa.o ~ 0.7) the 
critical shell currently has a dimensionless radius fo = 0.84, 
i.e. it has a value 84% of its maximum radius. 



2.3 Conditions for a critical shell 

For the translation of the radius fcs of a bound object into 
a density criterion, we evaluate the average mass density ps 
enclosed by a given shell, 

3M 



Ps = 



(12) 



47rr^ 

in terms of the critical density pc = SH^ /SnG. The shell is 
bound if 

2nA 



El > El 



2.36. 



(13) 

pc pc r% 

The value 2.36 corresponds to a Universe with Q,a,o ~ 0.7. 
Note that for JIa = 1 {t —> oo) the critical shell's density is 

pca/Pc = 2. 

The corresponding density excess 5 of the spherical 
mass concentration with respect to the global cosmic back- 
ground pu(t) (= pc,o^m,o/ci^ , with pc,o the critical density 
at the present epoch), may be inferred from the expression 



1 + 5 



2nA0 /aV 

nr. 



vo /ay 
,0 \rJ 



(14) 



On the basis of this equation we find that a critically bound 
shell at the present epoch has a density excess in the order 
of Scs ~ 6.9. 



2.4 Bound Object: 

linearly extrapolated density threshold 

To be able to identify the primordial regions that corre- 
spond to bound, collapsing and/or virialized objects at any 
arbitrary redshift z, we need the value of the corresponding 

linearly extrapolated densities. 

According to gravitational instability theory (jPeebled 

Il980ll . in the linear regime the density excess 5{a) grows as 



5{a) = D{a)So 



(15) 



where D{a) is the linear density growth factor (growing 
mode). In a FRW Universe with matter and a cosmologi- 
cal co nstant, it can be computed from (|Heath,1977. : .Peebles! 
Il980ll 



D(a) = ag{a) = H 



(a) r 

Jo 



da' 



(16) 



where g{a) is the growth with respect to that in an Einstein- 
de Sitter Universe {D{a) is normalized such that D{a) « a 
for a 0). 

In order to find the linear density excess So for the critically 
bound shell, one should evaluate its evolution at early epochs 
{a <^ 1). At these early times - when density perturbations 
are still very small, 5^1- the linearly extrapolated density 
excess (Eqn. llSp represents a good approximation for the 
(real) density of the object (Eqn. [13]). 

Using the fact that the early Universe is very close 
to an Einstein-de Sitter Universe and expands accordingly, 
a{t) cx t'^^^, we may infer that the early density excess 5 
of the bound sphere and its dimensionless radius r is (see 
appendix 



m = ^ m ■ 



(17) 



Using the approximate evolution of f (a) oc a implied by 
Eqn. (|A3|I . for early times (l-|-5«l)we find 



= To I a;:^; 



(18) 



The above leads us directly to the value of the linear density 
excess So, as at the early Einstein-de Sitter phase S{t) ~ 
a{t)5o, 



S = i ^ 1- /^ 2»A.o 

a 10 V flm.Q 



1/3 



1.504 . 



(19) 



The corresponding linearly extrapolated density excess for 
marginally bound structures at the present epoch is 

5t{a=l) = 1.17, (20) 

where we have used the approximation for g{a) = D{a)/a 



Future Evolution of Bound Superclusters in an Accelerating Universe 5 



l|Carroll et al.lll992l ) 

5 

g{a) « -n,„(a)x 



These objects are due to reach turnaround at a ^ oo. 



(21) 



2.5 Tests of the Spherical Binding Criterion 

Reahty is always far more complex than a simple an- 
alytical criterion is liable to cover. In addition to dis- 
tinct anisotropics, internal inhomogeneities and velocity dis- 
persions will influence the viability of the derived spher- 
ical binding c riterio n. By means of N-body simulations 
iDiinner et al.l (|2006l ) tested the binding density criterion 
Pcs/pc = 2.36 (for the current epoch) and the criterion in- 
volving the mass enclosed within the radius res. On average 
72% of the mass enclosed within res is indeed gravitationally 
bound to the structure. At the same time it was found that 
a mere 0.3% of the mass bound to the object is not enclosed 
within this radius. 



3 THE COMPUTER SIMULATION 

We simulate a standard flat ACDM Universe with cosmolog- 
ical parameters fim.o = 0.3, f^A.o = 0.7, and h — 0.7, where 
the Hubble parameter is given by Ho = lOO/i km s"^Mpc"\ 
The normalization of the power spectrum is as = 1. In or- 
der to have a large sample of bound objects, the simulation 
box has a side length of 500/i~^ Mpc and contains 512'^ 
dark matter particles of mass nidm = 7.75 X lO^Oft-'Mo 
(see Fig. [1} . The initial conditions are generated at expan- 
sion factor a — 0.02 (redshift z = 49), and evolved until 
a — 100 usin g the ma,ssive parallel tree N-Body/SPH code 
GADGET-2 (|Springelll20"05l ). The Plummer-equivalent soft- 
ening was set at ep; = 20 /i~^kpc in physical units from 
a = l/3toa = 100, while it was taken to be fixed in co- 
moving units at higher redshift. Given the mass resolution 
and the size of the box, our simulation allows us to reli- 
ably identify massive superclusters with ~ 80, 000 particles. 
The simulation was performed on the Beowulf Cluster at the 
University of Groningen. 

We took snapshots at five different timesteps: starting 
at the present time (a = 1), we studied the mass distribution 
at a = 2, a = 5, a = 10, to ultimately end up in the far 
future at a = 100. The latter was taken as a representative 
epoch at which the internal evolution of all bound objects 
appeared to have been completed. 



3.1 Identification of bound objects and 
superclusters 

We identify and extract groups and objects in the simu- 
lation on the basis of a random subsample of 256'^ parti- 
cles (1/8 of the total particle number). For the identiflca- 
tion of bound structures we apply a four-step procedure. In 
the first step, we find all virialized halos in the simulations 
with more than 50 particles. Subsequently, we incorporate 



the surrounding spherical region bound to these halos. We 
then join the bound spheres that overlap each other into 
single bound objects. Finally, among these objects we select 
the most massive ones, the superclusters in the simulation 
volume. 



3.1.1 Mass range 

A first practical issue is the minimal amount of particles we 
deem necessary to accept a halo/object detection as signif- 
icant. We choose a minimum of 50 particles, corresponding 
to a mass cut of M ^ 3.1 x IO^^/i'^Mq. This means our 
cluster and bound object sample has an implicit bias in not 
containing any objects with a mass less than the mass limit. 
This is not a problem for the most massive objects, but may 
produce an incomplete sample for lower mass objects. 



3.1.2 HOP and virialized halo finding 

In order to find gr oups of particles presen t in our simula- 
tion, we use HOP l|Eisenstein fc HutlflQQSl ). This algorithm 
first assigns a density estimate at every particle position by 
smoothing the density field with an SPH-like kernel using 
the Udens nearest neighbors of a given particle. In our case, 
we use Udens = 64. Subsequently, particles are linked by as- 
sociating each particle to the densest particle from the list of 
its Uhop closest neighbors. We use Uhop ~ 16. The process is 
repeated until it reaches the particle that is its own densest 
neighbor. 

The algorithm associates all particles to their local max- 
ima. This procedure often causes groups to fragment. To cor- 
rect this, groups are merged if the bridge between them ex- 
ceeds some chosen density thresholds. Three density thresh- 
olds are defined as follows (|Cohn et al.ll2001^ : 

- Sout'. 

the required density for a particle to be in a group. 

^saddle • 

the minimum boundary density between two groups 
so that they may have merged. 

the minimum central density for a group to be 
independently viable. 

We follow the criterion of lEisenstein fc Hull (|l998h : 

5outer.5saddie.3peak='i-.'2.5:3. The valuc of Speak is associatcd 
with that of the corresponding density Avir{a.) of the virial- 
ized core of the bound regions. 



3.1.3 Virial Density Value 

To determine the value of the virial density A„ir(a) in the 
HOP formalism we resort to the dynamical evolution of 
spherical top-hat perturbation. The value of A^ir is ob- 
tained from the solution to the collapse of a spherical top- 
hat perturbation under the assumption that the object has 
just virialized. Its value is Avir = 18n^ for an Einstein- 
de Sitter Universe. For the cosmology described here, at 
a = 1 its value is Avir{a — 1) x 337. The value at 
later epochs, or in other cosmologies, is obtai ned by solving 
the spherical collapse equations numericallv (iGunn fc Gott 
ll972l : lLacev fc Coldl 19931 : lEke et al.|[l99^ : iKitavama fc Suto 
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Figure 1. Simulated mass distribution at a = 1 (left column) and a = 100 (right column). The top panels show the particle distribution 
in a 500/i~^Mpc box, in a 30/i~^Mpc thick slice projected along the z-axis. Top left: mass distribution at a = 1. Top right: mass 
distribution at a = 100. Lower left: zoom-in on the a = 1 mass distribution within 50/i~^Mpc box indicated by the square in the top 
lefthand panel (similar 30/i~^Mpc thick slice along z direction). Lower right: region of same 50?t~^Mpc physical size as the lower left 
panel, at a = 100, around the same bound supercluster. 



1 19961 : iBrvan fc Norman 1998t). An exte nsive description of 
this can be found in lArava-Meld l|2008l '). 



At a = 100, Q,m = 4.3 x 10~^. For this situation we 
resort directly to the virial theorem to determine the char- 
acteristic virial radius. According to the virial theorem the 
kinetic energy K^ir of a body whose pot ential is of the form 
Vvir = R" is equal to K yjr = (n/2)Kir l|Landau fc Lifshit3 
Lahav et al.lll99ll ). For a general case of a virialized 
object in a Universe with matter and a cosmological con- 



stant, 



Kyir — ,vir ^A,vir 



(22) 



Note that here we write energies in dimensionless form (see 
e.g. Eqn. [5]and Eqn. [Sjl. In this equation Vg is the gravita- 
tional potential energy and Va is the potentia l energy due 
to the cosmological constant l|Lahav et aLlll99ll ). Hence, the 
total energy Evir ~ K^ir + Vvir of a virialized object is 



(23) 



Future Evolution of Bound Suverclusters in an Acceleratino Universe 7 



Because of energy conservation, the energy at maximum ex- 
pansion is equal to the energy at virialization. This trans- 
lates into the following equation for the dimensionless radius 



This is a cubic equation with solutions f ~ —1.366, f = 
1 and f ~ 0.366. The first is an unphysical solution. The 
solution f = 1 is the maximum radius of the critically bound 
shell (see Eqn.[23J. The third value, f « 0.366, is the virial 
radius of the enclosed mass M. The corresponding virial 
density excess of the mass clump would be 1 + Avir = p/Pc ~ 
40.8. 



3.1.4 Clusters 

Besides the definition of a sample of superclusters, we also 
identify the clusters in our simulation. Their identity is more 
straightforward to define, since we may presume they are 
virialized. 

Virialized halos, as identified by HOP (sect.[3TT2J, with 
masses larger than 3 x lO^'^/i^^M©, are considered clusters. 
Note that, by definition, these clusters are identical to the 
bound object cores produced in step 1 of our supercluster 
finding procedure (sect. [37f .2p . 

3.1.5 Bound Halo Identification 

Once we have identified the virialized halos in our sample, 
we have to proceed by outlining the gravitationally bound 
region around these cores and, if necessary, join them to- 
gether into a supercluster. In practice, we do this as follows. 

Of an identified IfOP halo, we take the location of the 
densest particle as a first estimate of the center of mass. 
Subsequently, we grow a sphere around this center, with the 
radius being increased until the mean overdensity within the 
corresponding radius reaches a value of 300pc. This value is 
chosen in order to find the densest core of the structure. We 
then calculate the center of mass of this sphere and repeat 
the process, iterating until the shift in the center between 
successive iterations is less than 1% of the radius. 

With the final center of mass, we apply the criterion of 
Eqn. [13] for identifying the bound spherical region around 
the HOP core. To this end we determine the radius rca at 
which the average interior density pes/ pc reaches a value: 
p^^/pe = 2.36 at a = 1, 2.22 at a = 2, 2.06 at a = 5, 2.02 at 
a = 10 and 2.00 at a = 100 (see sect. [23]). 

3.1.6 Joining Halos: Bound Objects 

The procedure outlined in the previous paragraph will fre- 
quently lead to overlapping bound spheres that in reality 
will be bound to each other. In order to account for this, 
we follow a radical prescription. If two spheres overlap, we 
proceed with the most massive one and join the lower mass 
sphere to the high mass one while removing it from the list 
of objects. 

We found that by following this procedure at a = 100, 
such overlaps do not occur. It implies that in the far future 
nearly all bound objects are compact and isolated islands 
in the Universe. At much earlier epochs, and specifically at 
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Figure 2. Scatter plots of the mass of identified bound objects, 
Mboundi against the mass of the virialized HOP halo that forms 
its core, Mhop- The density of points (MjjoPiMBound) in the 
scatter diagram is depicted in terms of an isodensity contour map. 
The highest concentration of points has the darkest (black) shade, 
gradually fading towards light colours. Top: a = 1. Bottom: a = 
100 

a = 1, it does turn out to occur for a significant fraction of 
the bound spheres. 

3.1.7 Superclusters and Bound Objects 

Because not all bound objects would be prominent enough 
to be a supercluster, we use a mass threshold to select the 
superclusters amongst the bound objects in our sample. The 
mass threshold is chosen to be Msc ~ 1O^^/!,~^M0, approx- 
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imately the mass of the Local Supercluster. Although this 
value is somewhat arbitrary, and a somewhat different value 
might also have been viable, it represents a reasonable order 
of magnitude estimate. 

As a result, we reserve the name bound objects for all the 
objects that have ended up in our sample, while superclusters 
are the subset with masses higher than our supercluster mass 
threshold of IO^'^Mq. 

3.1.8 Sample completeness 

Fig. [2] helps us to evaluate the completeness of our object 
sample. The panels contain scatter plots of the total mass 
of the bound objects versus the HOP mass of their corre- 
sponding virialized cores. 

Evidently, the lower right region is empty: HOP groups 
will always be less massive than the bound groups. There is 
a correlation between both masses, but with a high scatter. 
As expected, the evolution of bound objects towards fully 
virialized clumps expresses itself in a substantially stronger 
correlation at a = 100 than at a = 1. On the basis of these 
relations, we may conclude that at a = 1 the sample is 
complete for masses greater than 2 x 10^ h~^MQ, while at 
a = 100 the sample is complete for masses down to 6 x 
10^^/i~^Mq. In order to keep the samples comparable, we 
use a mass completeness threshold of 2 x 10^ /I'^M© at 
a = 1 and 1.4 x 10^'^h~^MQ at a = 100 (see discussion 
sect. E}. 

3.1.9 The Object Sample 

At a = 1 HOP finds ~ 20600 virialized "clusters" with more 
than 50 particles, i.e. halos with a total mass M ^ 3.1 x 
10^^/i"^Mq. At a = 100 it finds ~ 18000 virialized objects. 
These will be taken as the starting point for our supercluster 
finding procedure. They also constitute the cluster sample 
in our simulation. 

After determining the connected bound region and join- 
ing these overlapping bound spheres (see sec. I3.1.6p . plus 
checking for sample completeness, we finally end up with a 
sample of ^ 4900 bound objects at the current epoch. At 
the other epochs, from a = 2 to a = 100, this is approxi- 
mately the same number. Of these, ~ 535 are superclusters 
(at a — 1), while seventeen are truely massive supercluster 
complexes with Msc > 5 x lO^'^/i'^M©. 



4 EVOLUTION OF BOUND OBJECTS 

The two top panels of Fig.[T]show a slice of 30/i~^Mpc width 
of the particle distribution projected along the z axis, at 
a = 1 and a = 100. By taking a region of the same physical 
size at both epochs, the effect of the de Sitter expansion of 
the Universe becomes manifestly clear. 

At a = 1, the large-scale structure of the cosmic web is 
well established and its morphology and character hardly 
change thereafter. The lower left panel zooms in on the 
square region of the top-left panel. Centered on a massive 
structure, it shows the mass distribution at a = 1. The ra- 
dius of the circle is that of the bound region, according to 
the criterion of Eqn. 1131 It shows that it is well connected 
with the surrounding structure. The same object, but now 



at a = 100, is depicted in the lower right panel (with the 
same physical scale). We see that the size of the bound ob- 
ject is nearly the same at both expansion factors. While 
in comoving coordinates the accelerated expansion of the 
Universe results in a freezing of structure growth on scales 
much larger than the initial size of superclusters, in phys- 
ical coordinates the separation of structures continues and 
grows exponentially in time. This results in clearly detached 
regions which evolve in complete isolation: genuine cosmic 
islands. 



4.1 Case Studies 

To get an impression of the internal and external evolu- 
tion of the bound objects in our simulation, we focus on a 
few specific objects. By following qualitatively the "inter- 
nal" evolution of three bound objects we intend to set the 
scene for the further analysis in this study. 

Fig. |3] displays the evolution of one of the most massive 
objects in our sample. At a = 1, it has a substantial degree of 
substructure (Fig. 3] lefthand). The centre of the superclus- 
ter is dominated by the massive central cluster that is the 
virialized object from which we constructed the remainder of 
the supercluster. It forms the centre of a huge complex, con- 
necting the surrounding matter distribution via prominent 
filamentary extensions. These form the transport channels 
along which mass flows into the central supercluster region. 
Noteworthy is the large number and variety of subclumps 
along the fllaments and around the centre of the superclus- 
ter. This high mass supercluster undergoes a radical change 
towards the future. 

By a = 2 we see that most of the surrounding material 
has fallen into the central core, with a dramatic decrease 
in the number of surrounding subclumps. While at a = 5 
we still see a significant number of small clumps around the 
central supercluster, by a = 10 only a few individuals seem 
to have survived. Comparing the 1st and 2nd panel with 
the panels corresponding to a = 5 and a = 10 also seems to 
suggest that the infalling subclumps at later epochs have a 
lower mass. 

With respect to the morphological character of the sur- 
rounding mass distribution, we note that at a = 2 we can 
still discern the vague remnants of the salient filamentary 
patterns at the present epoch. Nonetheless, most of the for- 
merly richly patterned weblike structure has disappeared 
and seems to have resolved itself as its mass accreted onto 
the supercluster. The supercluster has also assumed a more 
smooth and roundish appearance, even though at a = 2 we 
can still recognize the original geometry, both in terms of 
its elongated shape and its orientation along the same di- 
rection. At even later epochs, the trend towards a highly 
centralized and regular mass concentration with a nearly 
perfect spherical shape continues inexorably. At a = 10 we 
may still just recognize some faint matter extensions at the 
edge of the supercluster. However, these do no longer bear 
the mark of its original orientation and shape. At the final 
timestep, a = 100, the object has reached the ultimate con- 
figuration of a perfectly regular, nearly spherical, and cen- 
trally concentrated and largely virialized dark matter halo. 
All surrounding substructure within the binding radius has 
fallen and has been entirely absorbed while the supercluster 
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Figure 3. Evolution of one of tlic most massive objects in our simulation, witii a present-day mass of 6.8 X lO^^'h^^MQ, at a = 
The supercluster is shown at a = 1, a = 2, a = 5, o = 10 and a = 100 (zigzagging from top left to bottom left). The size of the box 
always 14/i~iMpc in physical coordinates. 



10 Araya-Melo et al. 




Figure 4. Evolution of two different bound objects in the simulation. Each row shows the evolution of a single object at a = 1, a = 2, 
a = 5, a = 10 and a = 100. The particle positions are in physical coordinates. The box sizes have been scaled to the size of the mass 
concentrations and therefore differ amongst each other. Top row: intermediate mass bound object, M~ 5.6 X 10^''/i~^Mq (a = 1). Boxsize 
is 6/i~^Mpc . Bottom row: least massive bound object, with a present-day mass of ~ 2 X lO^*h~^M0 (boxsize at a = 1: 4/i~^Mpc, 
boxsize at a = 100: 3/i~^Mpc). The evolution of these medium and low mass objects should be compared to that of the the massive 
supercluster in Fig. |3] 



attained a perfectly virialized configuration: the big mass 
concentration has become a true island universe. 



THE SPATIAL DISTRIBUTION OF BOUND 
OBJECTS AND SUPERCLUSTERS 



The evolution of the massive supercluster is compared 
with that of two more moderate bound mass clumps in 
Fig. |4l The top row shows the evolution of a medium mass 
bound object, with a mass of M ~ 5.6 X lO"/i-^M0 (at 
a = 1), while the bottom row depicts the development of 
the least massive bound object. The present-day mass of 
the latter is M ~ 2 x lO^''/i~^M0. Even though the medium 
mass object (central row) shares a similar trend towards a 
centrally concentrated virialized clump, we also notice that 
its influence over the surroundings is considerably less pro- 
nounced and extends over a considerably smaller region. At 
o = 1 we recognize some relatively large subclumps in its 
surroundings, most of which by a = 2 have fallen in. The 
subclumps and the more diffusely distributed surrounding 
matter do not seem to display a pronounced spatial pattern 
and they do not appear to be organized in weblike filamen- 
tary extensions. Even though to some extent this may be 
a consequence of the limited resolution of our simulation, 
it undoubtedly pertains also to the substantially lower dy- 
namical (tidal) influence o f the clump over its surroundings 
(see e.g. iBond et al.1 ligQel ). After a « 2 - 3, the infaU of 
matter proceeds mostly through quiescent accretion, result- 
ing in a gradual contraction of the object into a moderately 
elongated ellipsoidal halo. 

Even less outstanding is the evolution of the low-mass 
bound object (bottom row). It hardly shows any substruc- 
ture and seems to consist only of a central region and a few 
particles within the binding radius. Between an expansion 
factor of a = 2 and a = 100 the changes in appearance 
are only marginal in comparison to those seen in the more 
massive bound structures. 



The evolving spatial distribution of the bound objects in 
our simulation is shown in Fig. [S] From this direct visual 
inspection, we see that the spatial distribution of the bound 
structures hardly changes between a = 1 and a — 100. 

Immediately striking is the fact that the bound object 
or supercluster "flUing factor" f dramatically declines as 
time proceeds. Defining the filling factor as the fractional 
volume occupied by bound objects and superclusters. 



f = ^Vj/Vt 



(25) 



we may appreciate the strength of the effect by inspecting 
Table[T] At the onset, ie. the present time, the bound objects 
and superclusters take up a sizeable fraction of the cosmic 
volume: ~ 3.3% for the bound objects and ~ 1.2% in the 
case of the superclusters. Even while towards later times 
the number of bound objects with viralized cores is increas- 
ing (unlike that for the superclusters), the volume fraction is 
rapidly declining. After a = 10 they occupy a negligible frac- 
tion of the Universe, superclusters only a fractional volume 
in the order of ~ 10~*. The tremendous exponential expan- 
sion of the Universe as a result of the cosmological constant 
clearly renders each bound cluster and supercluster an ever 
more isolated and lonely island in the Universe ! 

The other important observation is that the spatial dis- 
tribution of bound clumps hardly changes from a = 1 to 
a — 100. The spatial patterns visible in the object distribu- 
tion at a = 1 are very similar to the ones at a much later 
epoch. Clustering does not seem to get any more pronounced 
after the present epoch a = 1. 
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Figure 5. Spatial distribution of bound structures at five different expansion factors: a = 1 (top left panel), a = 2 (top center panel), 
a = 5 (top right panel), a = 10 (bottom left panel) and a = 100 (bottom right panel). The panels shows the distribution in comoving 
coordinates in a 100/i^~Mpc thick slice projected along the z-axis. The bound structures are presented as circles with a size equal to 
their comoving radius (note that the comoving size of the bound objects shrinks as the Universe evolves into the future). Amongst the 
bound objects, the dark circles indicate the superclusters (M^ IO^^/i^^Mq). 



5.1 Correlation Analysis 

To quantify the visual impression of the spatial distribution 
of bound objects and superclusters in our sample, we have 
determi ned their spat ial two-point correlation function ^(r) 
(see e.g. |PeeblesllT980l ) at the different epochs. 

The left frame of Fig. |6] presents the two-point correla- 
tion function ^(r) at a = 100 for three different samples: the 
complete simulation sample of bound objects, the subset of 
bound objects with mass M> lO^^M© and the superclusters. 
We assess ^(r) in terms of comoving distances, in physical 
coordinates one should evidently take into account the di- 
lution factor due to the expansion of the Universe. In our 
correlation analysis we take a = 100 as reference point, as 
we will see that the situation at other epochs is more or less 
comparable. In all three cases we observe the strong cluster- 
ing of the populations. Over a large range, almost extending 
out to a distance of r « 100/i~^Mpc, the bound objects and 
superclusters do show a significant level of clustering. This 
can be inferred from the inspection of the lin-lin plot of ^(r) 
in the insert in Fig. [S] (lefthand frame) . 

To first approximation, ^(r) behaves like a power-law 
function, 

= . (26) 



Table 1. Filling factor of bound objects (left) and superclusters 
(right) at 5 expansion factors: a=l, a=2, a=5, a=10 and a=100. 





Bound Structures 


Superclusters 


a = l 


3.34x10-2 


1.24x10-2 


a = 2 


4.34x10^3 


1.37x10-3 


a = 5 


3.31x10-"' 


8.64x10-5 


a = 10 


4.24x10-5 


1.01x10-5 


a = 100 


4.14x10-* 


0.91x10-** 



This is particularly evident for the sample of superclusters, 
which at a = 100 has a correlation length ro = 23±5/i-^Mpc 
and a slope 7 ~ 2.1±0.1. Within the error bounds the num- 
bers at the other timesteps are equivalent. While the overall 
behaviour of the correlation function of all bound objects, 
or that of the ones with M> IO^^Mq, does resemble that 
of a power-law we also find marked deviations. At the small 
distances where the clustering is strongest, r < S/i-'^Mpc, 
we find that ^(r) has a distinctly lower slope, 7 ~ 1.4 — 1.5. 
In the large-scale regime, i.e. beyond the clustering length 
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Figure 6. Two-point correlation function ^(r) of bound objects and of superclusters. The correlation function is shown in a log-log plot. 
Left: Two-point correlation function at a = 100 for three different samples: all bound objects in simulation sample (blue dots), bound 
objects with mass M> 10^'*Mq (green dots), superclusters (red stars). The red line is a power-law fit to ^(r) of the superclusters. The 
insert is the lin-lin plot of ^(r) for the whole object sample. Right: the evolution of i;{r). The two-point correlation function of the whole 
objects sample at five different epochs, a = 1, 2, 5, 10 and 100. The corresponding colours are indicated in the legenda the top-right 
corner of the frame. For guidance the red line shows the power-law fit to the supercluster ^{r) at a = 100. 



ro ~ 11. 5h ^Mpc, the slope steepens considerably and at- 
tains a value close to 7 ~ 2.0 at o = 100. 

We find that more massive objects are more strongly 
clustered, confirming the impression obtained from Fig. [S] 
(see e.g. top-left frame) . Comparison between the clustering 
of the whole sample with that of the subsample of objects 
with mass M> lO^^Mo, and in particular that of superclus- 
ters, shows that the more massive samples have a stronger 
correlation function over the whole range of distances. In 
addition to this, we may also observe that the clustering of 
superclusters seems to extend out to larger distances than 
that of all objects. That is, for superclusters ^(r) > out to 
distances larger than r « 100/i~^Mpc. 

This clustering scaling behaviour is a telling illustra- 
tion of the (clustering) bias of more massive objects with 
respect to the global average, and has to some extent be 
related to the properties of the peaks in the Gaussian pri- 
mordial density field (Kaiser 198J). Our findings relate di- 
rectly to the known clus tering behaviour of clusters. Since 
iBahcall fc Soneiral ()l988i ) found that the clustering strength 
of clusters is an increasing function of the cluster richness, 
a large variety of studies have uncovered the systematics 
of this clustering bias (see e.g. ISzalav fc Schramml 1 19851 : 
iPeacock fc West! il992l : lEinasto et all |2002|) . In particular 
enticing are the recent results obtain ed for cluster sam- 
ples extracted frorn the SDSS survey (|Bahcall et al.l [200^ : 
lEstrada et"aLll2009^ . 



In the righthand frame of Fig.[6]we assess the evolution 
of the two-point correlation function of the complete sample 
of bound objects, from a = 1 to a = 100. With the exception 
of some evolution at short distances, r < 5/i~^Mpc, there 
is not any significant evolution of the two-point correlation 
function over the full range of distances. This is a mani- 
fest illustration of the stagnation of structure evolution in 
an accelerating Universe on scales larger than several Mega- 
parsec. The situation appears to be exactly the same for the 
stronger correlation function of superclusters: from a = 1 to 
a — 100 superclusters keep the same level of clustering. Of 
course this is a clear manifestation of the expected end of 
structure growth after the Universe went into acceleration 
at 2; « 0.7. Note that for reasons of clarity, in the righthand 
frame of Fig. |6] we restrict ourselves to merely plotting the 
power-law fit to the supercluster correlation function. 



6 GLOBAL EVOLUTION: MASS FUNCTIONS 

With the identification of superclusters we may ask how 
much mass is contained in them at the present epoch and at 
a — 100. At a = 1, the entire sample of bound objects in our 
simulation box amounts to 2.73 x 10^**/i~^Mo. At a = 100 
this has grown to a total mass of 2.83 x 10^*/i~^Mq. This 
represents 26%, respectively 27%, of the total mass in our 
simulated Universe (1.04 x lO^'^/i'^Me). 
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We may make two immediate observations. First, given 
that bound objects tend to lose around 28% of their mass 
towards a = 100 (see below), it must mean that the number 
of bound objects fulfilling our criterion - of containing one 
or more virialized cores - is still growing from a = 1 towards 
a = 100. It may also indicate a problem in applying a purely 
spherical density criterion: at a = 1 mass concentrations are 
more aspherical and inhomogeneous, while at a = 100 they 
are nearly spherical concentrated mass concentrations (see 
sec. [7]). Perhaps even more tantalizing is the fact that ap- 
parently more than 70% of mass in the Universe will remain 
outside of the supercluster islands and will keep on floating 
as a lonely population of low-mass objects in a vast cosmic 
void. 

Perhaps the most outstanding repercussion of the slow- 
down of large-scale structure formation in hierarchical cos- 
mological scenarios is the fact that the condensaton of new 
objects out of the density fleld will gradually come to a halt. 
This should be reflected in the mass spectrum of the objects 
that were just on the verge of formation around the time 
of the cosmological transition. Here we investigate the mass 
distribution of superclusters, i.e. bound but not yet fully 
virialized structures. For comparison we also investigate a 
sample of virialized halos. 



6.1 Theoretical mass functions 

Fig. [7] shows the evolution of the mass function of viri- 
alized objects predicted by the Press-Schechter formalism 
l|Press fc Schechteilll974) and the Sheth-Tormen excursion 
set prescription ( Sheth fc Tormenlll999l ) , for the current cos- 
mology. The Sheth-Tormen expression takes into account 
the anisotropic collapse of dark halos. We also compare 
our mass functions t o the heuristic simulation-based mass 
function suggested bv I Jenkins et al.l l|200ll '). We refer to ap- 
pendix[B]for a listing of the expressions for these mass func- 
tions. 

The top panel shows the strong evolution in the past, 
by depicting the mass functions at 2 = 4, 2, 1, 0.5 and 0. It 
shows that structure grows in mass and number while the 
Universe expands. The future evolution is a lot less strong, 
as evidenced by the mass functions at expansion factors a — 
1,2,4,10 and 100 in the bottom panel. After a = 1 the 
number of low mass objects does not change substantially, 
while after a = 4 evolution comes to a complete halt. This 
may be best appreciated from the fact that the curves for 
a = 10 and a = 100 overlap completely. 

We see that the Sheth-Tormen approximation predicts a 
higher number of massive clusters than the Press-Schechter 
formalism. As anisotropic collapse speeds up the contraction 
along the minor axis of an object, there is a higher number of 
regions reaching a sufficiently large overdensity before dark 
energy prevents any further evolution. Implicitly, this lowers 
the number of low-mass objects as more get absorbed into 
the high mass superclusters. 



6.2 Simulation mass functions 

We assess the mass function of the virialized halos, the ob- 
jects identified by HOP on the basis of the prescription in 
sec. 13.1.21 and that of the bound (supercluster) objects that 
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Figure 7. Theoretical mass functions for virialized objects in a 
ACDM Universe (r2m,o = 0.3, Oa.o = 0.7, /x = 0.7). Shown are 
the Press-Schechter mass functions (dashed lines) and the Sheth- 
Tormen mass functions (solid lines). Top: the mass functions at 
redshifts z = 4,2,1,0.5 and (from left to right). Bottom: the 
mass functions at expansion factors a = 1, 2, 4, 10 and 100 (from 
left to right). 



were identified according to the criterion specified in sec- 
tion [3TT] Given that the we do not expect a radical change 
in mass functions between the present epoch and a — 100, 
we restrict the comparison to the simulation mass functions 
at a = 1 and a = 100. 

Fig. |8] shows the mass functions of the virialized halos 
found by HOP (left panel) and of the bound objects (central 
panel), at a = 1 and a = 100. As expected, the number of 
massive virialized halos increases as we go from a = 1 to 
a = 100. The increase is only minor, yet significant, and 
a manifestation of the freezing of structure formation and 
hence that of the corre sponding mass functions (see also 
iNagamine fc Loebllioosl ). The mass functions at a = 1 and 
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Figure 8. Mass functions of virialized halos and of bound objects. Left: integrated mass function N(>M) for virialized HOP objects 
in the simulation, at a = 1 (solid line) and a = 100 (dotted line). Centre: integrated mass function N(>M) for bound objects in the 
simulation, at a = 1 (solid line) and a = 100 (dotted line). Right: the integrated mass function N(>M) for bound objects compared with 
the reduced a = 1 mass function (dashed line). See text for explanation. 
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Figure 9. Mass functions of virialized halos and of bound objects compared with three theoretical mass functions. These are the Press- 
Schechter mass function (dotted line), the Sheth-Tormen function (dashed line) and the Jenkins function (dot-dashed line). For the 
virialized HOP halos the critical overdensity 5c for the PS and ST functions is the one for collapse, for the bound objects the value 5^, for 
assuring a bound object (Diinner et al. 2006, this study). Left panel: the a = 1 integrated mass function N(>M) of HOP halos. Centre: 
the a = 100 integrated mass function N(>M) of HOP halos. Righthand: the a = 100 integrated mass function N(>M) of bound objects. 



a — 100 also reflect the continuing hierarchical evolution 
within the realm of the bound supercluster regions. There 
is a definitive increase in the number of the most massive 
clumps, going along with a decrease at the low-mass side of 
the mass function. 

When turning towards the mass function for the bound 
(supercluster) regions (central panel) we flnd the same mass 
function at a = 1 and a — 100, except for a slight decrease 
in the mass of the objects over the whole mass range. There 
is a loss of mass, amounting to some 28 ± 13% of the mass 
enclo sed within the critic al radius of the superclusters at 
a = 1 iDiinner et al.l l|2006l ). This is substantially more than 
the mere 1% mass gain as a result of accretion of mass in 
between a = 1 and a = 100. The loss of mass has to be as- 
cribed to the virialization process of and within the bound 
object. A major factor in this is the abundant substructure 
in the supercluster at the present epoch as opposed to the 
smoothened mass distribution within the ultimate super- 



cluster island at a = 100. In order to correct for this "re- 
duced" mass we renormalize the supercluster mass function 
into a reduced supercluster mass function, simply by multi- 
plying the masses by a factor 0.72. As may be observed in 
the righthand panel of Fig. [8] we find an almost perfect over- 
lap between the mass function at a = 100 and its reduced 
equivalent at a = 1. 



6.3 Comparison of simulated and theoretical mass 
functions 

Fig. [9] shows the cumulative mass function of the virial- 
ized objects found by HOP at a = 1 and at a = 100, and 
compares them with three theoretical mass functions (see 
app.m. 



Future Evolution of Bound Superclusters in an Accelerating Universe 15 



6. 3. 1 Mass functions of hound objects and superclusters 

To compare the mass function of the bound (superclus- 
ter) regions with the theoretical Press-Schechter and Sheth- 
Tormen functions we need to specify a critical overdensity 
corresponding to the bound (supercluster) regions in our 
sample (see app.|Bj. Because in this study these are assem- 
bled on the presumption that they are marginally bound, we 
use the corresponding value of the linear extrapolated den- 
sity excess, 5h = 1.17, derived in section [23] (see Egn. I20p . 

At a = 1, the Sheth-Tormen function seems to provide 
a better fit than the Press-Schechter function to the bound 
object mass spectrum, in particular for the tail of massive 
superclusters. It demonstrates the importance of morpholog- 
ical and tidal influences on the mass spectrum of these gener- 
ically nonspherical objects (see next section [7]). At a = 100 
the Press-Schechter function provides a substantially better 
fit to the bound object mass function. This may be related 
to the fact that the PS formalism implicitly assumes pure 
spherical collapse, which we will see in sect. [7] agrees quite 
well with the shape of superclusters at a = 100. Also, we find 
that in the far future, when most of our superclusters have 
become isolated and largely virialized islands, it is appro- 
priate to compare with a Press-Schechter or Sheth-Tormen 
mass function based on another critical density value. In be- 
tween a = 1 and a = 100 most cosmic islands will reside in a 
dynamical phase somewhere between marginally bound and 
full collapse, implying a critical value in between 5b — 1.17 
(Eqn.[2Sl) and the collapse threshold Sc = 1.675 (Eqn.|B4|. 

It is also clear that the Jenkins function does not pro- 
vide a suitable fit to the supercluster mass function. This 
may not be surprising given the fact that it is a numerical 
approximation of the mass function of collapsed and viri- 
alized halos in N-body simulations and as such does not 
explicitly include an adjustable density threshold Sc. 



Neither the pure Press-Schechter nor the Sheth-Tormen 
function manage to fit the mass spectrum at a = 100 over 
the entire mass range. Press-Schechter does agree at the high 
mass end while Sheth-Tormen results in a better agreement 
at lower masses. This may be an indication for the more sub- 
stantial role of external tidal forces on the evolution of the 
low mass halos. Such external influences are entirely ignored 
by the Press-Schechter formalism, while they are succesfuU y 
modeUed by the Sheth-Tormen fits (|Sheth fc Tormenlll999l ). 



7 SHAPES OF BOUND STRUCTURES 

While the formation and evolution of structure on large 
scales comes to a halt once the Universe starts to accel- 
erate, the internal evolution of overdense patches continues. 
One of the most telling manifestations of the internal evolu- 
tion of these collapsing objects is their changing shape. The 
substructures that are within the bound radius will merge 
with each other into an increasingly smooth and concen- 
trated clump that will gradually assume a more and more 
spherical configuration. 

Using a variety of definitions for superclusters, 
their shape has b e en st ud ied both using real data 



Plionis et al] 1 19921: ISathvaprakash et al.1 Il998l: 



Einasto et al.l 



Basilakos et al. I I2OOII : ISheth. Sahni fc ShandarinI |2003| : 



■ .2007(3) a nd in N-Body simu l ations 

(e.g. ^Sathvaprakash ct al.' '1998'; Shandarin et al. '200^; 
Basilak os ot al. 2006; Wray ct al, 2006 ; Einasto e t al] 



Most studies agree that the dominant shape of 
superclusters at the present time is prolate, most evident in 
the presence of elongated filaments. These predominantly 
anisotropic shapes are a clear indication for the quasi- 
linear dynamical stage at which we find the present-day 
superclusters. 



6.3.2 Cluster mass function 

It is more straightforward to compare the mass function of 
the clusters, or HOP halos, in our simulation with that of 
the three theoretical mass functions for virialized objects. 
For these fits we use the critical collapse density value Sc = 
1.675, the value for our cosmology according to Eqn. IB4I 
The Jenkins approximation is of course independent of the 
value of Sc. 

At a = 1 the Jenkins mass function (Ijenkins et al.|[200ll ) 
is the one that fits best (dot-dashed line), perhaps not en- 
tirely surprising given its N-body simulation background. 
The Press-Schechter mass function represents a good fit at 
the lower mass end, al though it underestimate s the number 
of high mass clusters. iGovernato et al.] l|l999l ) claim that a 
critical density value Sc — 1.775 would provide a better fit 
and indeed it would lead to a small yet significant improve- 
ment of the Sheth-Tormen mass function. 

For a = 100 we adjust the parameters of the Jenkins 
function, using the fitting parameters for i}m = listed in 
lEvrard et al.l (|2002 i). With these parameters, it agrees very 
well with the HOP mass function, although it slightly over- 
estimates the number of lower mass objects. However, for 
Jenkins' original parameter values, it would lead to a sig- 
nificant overabundance of objects with respect to the ones 
found in the simulations. 



7.1 Definitions 

In order to determine the shape, we calculate the inertia ten- 
sor using all particles inside the spheres defined by Egn llSI 
with respect to their center of mass: 

lij = ^ XiXjm . (27) 

Since the matrix is symmetric, it is possible to find a coordi- 
nate system such that it is diagonal, yielding the eigenvalues 
ai, 02 and 03. These give a quantitative measure of the de- 
gree of symmetry of the distribution. With major axis si, 
medium axis S2 and minor axis S3, ie. S3 < S2 < si, the two 
axis ratios S2/S1 and S3/S1 are given by 




where > a2 > a-^. 

The object is almost spherical if both ratios S2/S1 and 
S3/S1 are close to one. Oblate objects have axis ratios S3 <C 
S2 ~ si, prolate objects S3 ~ S2 ^ si. 

7.2 Shape evolution 

Fig. \W\ reveals the overall evolving shape distribution in 
from the present epoch until a = 100. The shape distribu- 
tion is shown in terms of a plot of the axis ratios S2/S1 versus 
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Figure 10. Distribution of axis ratios for bound objects at five different expansion factors. Top row: a = 1, 2 and 5; bottom row: a = 10 
and a = 100. Each diagram shows the probability density distribution in the plane of S2/S1 versus ss/si values, with principal axes 
S3 < S2 < si. Plotted is the axis ratio probability distribution in iso-probability greyscales values, with percentile steps of 10%. Dark 
colors correspond to high probability, decreasing to lower probability as colours fade to light (in stops of 5%). 



Table 2. Average values of the axis ratios S2/s\ and sg/si, with 
their standard deviation, at five expansion expansion factor: a = 
1, a = 2, a = 5, a = 10 and a = 100. 





(S2/Sl> 


(ss/si) 


"■82/31 




a = 1 


0.69 


0.48 


0.13 


0.11 


a = 2 


0.73 


0.61 


0.14 


0.13 


a = 5 


0.87 


0.77 


0.10 


0.10 


a = 10 


0.90 


0.81 


0.08 


0.08 


a = 100 


0.94 


0.85 


0.03 


0.05 



ss/si of the objects. The corresponding shape distribution 
function is shown by means of greyscale iso-probability con- 
tour maps. Because by definition S3 < S2 < si, only the 
righthand triangle of the diagram is populated. Spherical 
groups are located at (1,1), oblate groups tend towards the 
fine S2/S1 = 1 while prolate groups are found near the diag- 
onal S2/S1 = S3/S1. 

At the first two timesteps, a = 1 and a = 2, we find 
a clear dominance of nonspherical objects. Most objects are 
distinctly anisotropic, which agrees with the observations 
of individual objects in sect. |4] At a = 1 the objects oc- 



cupy a wide range of mostly triaxial shapes, with a mean 
value ((S2/S1), {s3/si))= (0.69,0.48) and a standard devia- 
tion of ((T52/5J, (T3g/3j^)=(0.13,0.11) (see table[2ll. There are 
hardly any thin pancake-shaped structures, as low values of 
sa/si < 0.25 seem to absent. To some extent this may be a 
refiection of our bound group identification procedure, given 
its bias towards spherical configurations (see sect. I3.1.6p . 
Starting at a = 2 we observe a stretching of the shape dis- 
tribution (Fig. llOf) . along the direction of the prolate con- 
figurations, S3/S1 = S2/S1, and towards less anisotropic and 
more spherical morphologies. A sizeable fraction of the ob- 
jects shows a tendency towards a prolate shape. 

The fact that there are almost no spherical, or even 
nearly spherical objects, at the earlier epochs is hardly sur- 
prising given the relative youth of these supercluster objects. 
It is a refiection of the fact that t he primordial density field 
does not contain spherical peaks (|Bardeen et al.|[l986l ) and 
that the first stages of gravitational contraction and col- 
lapse proceed via strongly fiatt ened and elongated geome- 
tries (|Zeldovichlll970l : Ilck3ll973l ) . 

After a = 2 there is a rapid evolution of the vast ma- 
jority of objects towards a nearly spherical shape. At a = 5 
nearly all clumps are found in the corner with axis ratio 
S2/S1 = S3/S1 = 1. Also after a = 5 the width of the shape 
distribution continues to shrink. Not only do the bound 
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Figure 11. Distribution of axis ratios for three mass ranges. I: high-mass bound objects with M> 4.8 X lO^^/i ^Mq. II: medium mass 
range bound objects, 2.8xlO^'^h-^MQ < M ^ 4.8 X lO^/i-^M©. Ill: low mass bound objects, 2x lO^/i-^Mo < M ^ 2.8 X IO^^/i-^Mq . 
Each diagram shows the probability density distribution in the plane of S2/S1 versus sg/si values, with principal object axes S3 < S2 < «!■ 
Dark colors correspond to high probability, decreasing to lower probability as colours fade to light, with level steps corresponding to 
percentile differences of 10%. Top row: axis ratio distribution at the current epoch, a = 1. Middle row: axis ratio distribution at a = 5. 
Bottom row: axis ratio distribution at a = 100. 



clumps and superdusters on average attain a more spher- 
ical shape as the expansion of the Universe accelerates: it is 
the entire population of objects that appears to follow this 
trend. This systematic change of shape can be most clearly 
appreciated from inspection of table (2] 



a = 1 and the one at a = 100 reveals a manifest and even 
radical internal evolution of the supercluster complexes, at 
the same time when evolution on larger scales has been vir- 
tually frozen. 



In all, the contrast between the shape distribution at 
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Table 3. Number of objects and average values of the axis ratios S2/s\ and sa/si, with their standard deviation, at a = 1, a = 5 and 
a = 100 for three mass intervals (with mass range set at a = 1). 









(«2/si> 


{ss/si) 






I 


M > 4.8 X lO^-i/i-^Mo 


a = 1 


0.70 


0.49 


0.13 


0.11 






a = 2 


0.72 


0.60 


0.13 


0.12 






a = 5 


0.86 


0.76 


0.10 


0.10 






a = 10 


0.90 


0.81 


0.07 


0.08 






a = 100 


0.94 


0.85 


0.03 


0.05 


II 


2.8 X 10i*h-iMQ < M < 4.8 X lO^/i-iM© 


a = 1 


0.70 


0.48 


0.13 


0.11 






a = 2 


0.74 


0.63 


0.14 


0.13 






a = 5 


0.87 


0.77 


0.09 


0.10 






a = 10 


0.91 


0.81 


0.07 


0.08 






a = 100 


0.94 


0.85 


0.03 


0.05 


III 


2 X 1014/i-iMq < M < 2.8 X 10"/i-iMq 


a = 1 


0.68 


0.47 


0.14 


0.11 






a = 2 


0.73 


0.62 


0.15 


0.14 






a = 5 


0.86 


0.76 


0.11 


0.11 






a = 10 


0.90 


0.80 


0.08 


0.08 






a = 100 


0.93 


0.85 


0.04 


0.05 



T.3 Mass dependence 

One potentially relevant issue concerns the possible depen- 
dence of shape on the mass of bound structures. In or- 
der to investigate this, we divide our sample in three mass 
ranges, all approximately including the same number of 
bound structures. Class I are the most massive third of ob- 
jects, which at a = 1 have a mass M > 4.8 X 10"/j-^Mq. 
The medium mass class II includes objects with masses at 
a = 1 of 2.8 X lO^^h-'^MQ < M ^ 4.8 x 10"/i"^Mg and 
the low mass class III are those with masses at a = 1 of 
2 x lO"/i-^M0 < M ^ 2.8 X lO"/i-^M0. 

We found that over the entire interval of a = 1 to 
a = 100 there is hardly any distinction between the shape 
distribution in the different mass ranges (see Fig. Ilip . The 
three mass ranges have similar mean axis ratios at all five 
expansion factors (see Table |3}. This must be related to the 
fact that even though of different mass, the objects have 
been selected on the basis of similar (over)density values. 
The latter is an indication for the evolutionary state of the 
objects. In turn, as we have seen, this is reflected in their 
shape. 

We conclude that irrespective of their mass, all objects 
evolve into single, virialized and spherical objects as they 
grow in complete isolation after the Universe assumed an 
accelerated expansion and after all surrounding clumps and 
substructure within their binding radius has been accreted. 

7.3.1 Cluster vs. Supercluster Shapes 

Keeping in mind that the objects in our sample are bound 
but perhaps not virialized, it is instructive to contrast them 
to virialized galaxy clusters. Although their masses may be 
comparable, the bound group radii are much larger as they 
have substantially lower densities than clusters. Within this 
larger region there is a considerably more pronounced sub- 
structure. 

At a = 1, the clusters in our simulation are triaxial, 



with average axis ratios ({s2/si),{s3/si)) = (0.83,0.71) and 



standard deviation {a^ 



,3/,J= (0.09,0.09). These val- 



ues arc somewhat more pronounced than those quoted in 
other studies (e.g., Dubinski fc Carlbe r g 1991; Katz 199J; 
van Haarlem fc van de Weygaertl 1 19931: Ijir^^Suto 20o1: 



Kasun fc Evrardll2005l : [Paz et al.ll2006l : 



AUgood et al...2006i 'l. 



All agree that they tend to be more pr olate as the halo mass 
increases. iDubinski fc Carlberd l|l99lf ) found that halos are 
"strongly triaxial and very flat", with mean axis ratios 
of (s2/si)=0.71 and {s3/si)=0.50. For simula tions of iso - 
lated halos with different power spectra indices iKat j (|l99ll ) 
found values ranging from S2/si« 0.84 — 0.93 and sa/ai ~ 
0.43 — 0.71, while for massive clusters iKasun fc EvrardI 
(2005) found peak values of (s2/si, sa/si) = (0.76,0.64). 

Even though the radii of the bound structures in our 
sample are considerably larger than the virial radii of clus- 
ters, they do have similar triaxial shapes. Nonetheless, we 
do notice a distinct tendency of the superclusters to be more 
anisotropic than that of the virialized clusters. This must be 
a reflection of their different dynamical states. 



8 INTERNAL MASS DISTRIBUTION AND 
DENSITY PROFILES 

We have seen that superclusters are gradually decoupling 
from the global cosmic expansion, and towards the far fu- 
ture evolve in isolation as genuine cosmic islands. As these 
contract and collapse, and finally even virialize, they will 
develop into a much more compact object. To understand 
their internal evolution it is important to look into the devel- 
opment of their internal mass distribution, and hence their 
density profile. 
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Figure 12. Top: supercluster # 8 (M ~ 5.4 X lO^^h-^MQ). Bottom: supercluster # 98 (M ~ 3.6 X lO^^h'^Mo). Dark matter/particle 
distribution at a = 1 (left) and a = 100 (centre) of object # 8, in physical coordinates. For comparison the half mass radius (solid circle) 
and the virial radius (dashed-dot circle) are superimposed on the matter distribution. Righthand panel: Cumulative mass distribution of 
the object as a function of its (normalized) radius r/r;,, at a = 1 (solid line) and at a = 100 (dotted line). The cumulative mass M/M;, is 
normalized with respect to the corresponding final (bound) mass M;,. The solid vertical lines indicate the value of the half-mass radius 
(for which M/Mi,=0.5, the horizontal solid line), at both epochs a = 1 and a = 100. 



Figure 13. Radial density profiles of Object 8 (top panel) and Object 98 (bottom panel). In a log-log diagram each of the panels shows 
the radial density profile p(r), in units of the critical density Pcriti as a function of the normalized radius r/r^ir (normalized wrt. the 
supercluster virial radius). Solid line: profile at a = 1. Dotted line: profile at a = 100. For comparison in both panels we include a short 
line with an isothermal slope —2. 



8.1 Supercluster Case Studies 

Before trying to draw some general conclusions we look in 
more detail at the evolving mass distribution inside and 
around two representative individual superclusters in our 
sample. These are objects # 8 and # 98. The first one is 
a massive supercluster with a mass of ~ 5.4 X 10^^/i"^Mq, 
while the second one has a mass of ~ 3.6 X lO^^fe-^M©. Su- 
percluster # 8 is one of the most massive objects identified 
at a = 1 and also ends up as such at a = 100. 

The evolution of the mass distribution in and around 
superclusters 8 and 98 from a weblike irregular and struc- 
tured pattern at a = 1 into smooth and highly concentrated 
and nearly spherical mass clumps at a = 100 can be ob- 
served from the changing particle distributions in the left- 
hand and central frames of Fig. 1121 For guidance, on all four 
panels we superimposed circles centered at the superclus- 
ter's core. These are the "half-mass sphere" (solid circle), 
enclosing half of the total mass of the supercluster, and the 
"virial sphere", enclosing the central virialized core of the 
supercluster (dashed-dotted circle). 

For the structure of the virial core we turn to the log- 
log diagrams of the radial density profiles p{r) in Fig. 1131 
All profiles have the same basic shape, a high density cen- 
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Figure 14. Concentration parameter c of each bound object in the bound object sample as a function its mass M. The concentration 
parameter is the ratio of the half-mass radius of the object to its total binding radius rj, c = r^m/rj,. Superimposed on each of the plots 
is the line indicating the median in mass bins, along with the corresponding 25 and 75 percentile lines. Lefthand panel: a = 1. Righthand 
panel: a = 100. At a = 100 nearly all objects have a highly concentrated mass distribution, ie. c < 0.2. Note that the bound object 
sample is different at each expansion factor (see sec. I3.1.9t . 



tral core embedded within an isothermal power-law region 
with slope ~ —2. While there may still be some deviations 
from this slope at a = 1, comparison with the inserted short 
line of slope —2 shows that at a = 100 the object cores are 
almost perfectly isothermal. The profiles confirm the impres- 
sion from Fig. [T5] of the growth of the virial core, given the 
smooth near power-law profiles at a = 100 as opposed to 
the more irregular behaviour at the outer edges in the a = 1 
profiles. 

In addition, we note a radical change of the cosmic sur- 
roundings of both superclusters. At a = 1 both superclusters 
are still solidly integrated and embedded within the Mega- 
parsec Cosmic Web. Their central cores, indicated by the 
"virial spheres", are connected to the surroundings via fil- 
amentary tentacles along which we find a large variety of 
mass clumps. Note that these outer structures at a = 1 are 
actually bound to the supercluster core and is will fall in and 
merge with the central cluster as time proceeds. At a = 100 
the resulting supercluster concentrations have turned into 
isolated islands. 

It is telling that at a = 100 both superclusters have at- 
tained an almost equivalent internal mass distribution, even 
though at a = 1 their morphology was quite different. At the 
present epoch, supercluster # 8 is already a centrally con- 
centrated object. Given that the half-mass and virial spheres 
nearly overlap, we see that at a = 1 its virialized core con- 
tains nearly half of its total mass. This impression is under- 
lined by the cumulative mass distribution in the righthand 
panel of Fig. 1121 Comparison between the cumulative mass 
distribution at a = 1 (solid line) and a — 100 (dashed line) 



shows that there is only a small increase in mass in the inner 
region of the supercluster, immediately around the virialized 
core. Also, the fact that the half-mass radius moves only 
slightly inward implies a moderate change in the inner mass 
distribution. It is the outer half of the supercluster's mass 
which rearranges itself more strongly: as it falls in towards 
the supercluster's interior, the mass distribution becomes 
more concentrated and more regular. 

The changing mass distribution is more pronounced for 
supercluster # 98. At the current epoch its mass distribu- 
tion is more extended: its half-mass radius is located at an 
outward position. The central core is not nearly as promi- 
nent as that in supercluster # 8. The supercluster's mass 
increases rather slowly until the half-mass radius. Beyond 
this radius there is an abrupt rise until the outer superclus- 
ter radius. This is related to the presence of another major 
mass clump near the outer boundary. When we would have 
observed this supercluster in the observational reality, we 
would find it to be dominated by two very rich clusters. 

Also, we note that the superclusters do hardly gain mass 
from beyond their (binding) radius Rt- This is refiected in 
the flattening of the cumulative mass curves at a = 100. 



8.2 Supercluster Mass Concentration 

From the discussion in the previous subsection we have 
learned that the superclusters are turning into bodies with 
a highly concentrated mass distribution at a = 100. It is 
only towards these later cosmic epochs that the superclus- 
ters have turned into highly nonlinear and concentrated re- 
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gions. It would not be appropriate to seek to fit a theoreti- 
cal density profile to their radial mass distribution in order 
to determine their concentration. Instead of seeking to fit 
an Einasto profile (iEinasto ,1965i) , or the profusely pop ular 
universal NFW profile l|Navarro. Frenk fc White! Il997fl . we 



therefore prefer to define a concentration parameter that is 
independent of assumptions about the dynamical state of 
the mass concentration, 



(29) 



where rhm is the radius that encloses half of the mass. Note 
that with this definition a mass distribution with a "delta 
peak" in the centre would have c = 0, while a perfectly 
uniform distribution would have c ~ 0.8. An isothermal dis- 
tribution would correspond to c = 0.5. 

Fig. [14] shows the distribution of c as a function of mass, 
for the object population at a = 1 until a — 100. At the 
present epoch there still is a considerable spread of the con- 
centration parameter (c — 0.35, CTc — 0.14). There are even a 
few objects that get close to the c — 0.8 value corresponding 
to uniform mass profiles. Although there is a slight tendency 
towards higher concentrations, in general the concentration 
parameter at a = 1 refiects the irregular and prominent 
outer mass distribution. Low-mass bound objects appear to 
be more strongly concentrated than the A/ > 10^^ h~^MQ 
superclusters. 

The additional panels to reveal the expected develop- 
ment. After a = 1, there is a relatively rapid evolution 
towards much more concentrated configurations. While at 
a = 2 there are still a sizeable number of bound objects and 
superclusters with a rather diffuse or uniform internal mass 
distribution, at a = 5 nearly all objects have attained a high 
level of concentration of c < 0.4, with most having a concen- 
tration index c ~ 0.2. This development continues towards 
a radical conclusion at around a = 10. Over the whole mass 
range objects are highly concentrated from a = 10 onward. 
Also note that the trend of low-mass objects to be more con- 
centrated than high-mass superclusters persist, even though 
it gradually weakens towards later epochs. At the final time 
of a = 100, nearly all objects have c < 0.2, with an average 
concentration index c = 0.16 and (Tc = 0.02. 

It is clear that this evolution ties in strongly with the 
tendency to assume a near perfect spherical shape and, as 
we will see in the next section, with a rapidly decreasing 
level of substructure and multiplicity of the superclusters. 
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Figure 15. Supercluster multiplicity. The (dark) mater distribu- 
tion within the binding radius of a supercluster, at a = 1. The 
clusters within the supercluster realm are indicated by circles. 
The size of the circle reflects the (virial) size of each cluster halo. 



population within the supercluster is indicated by circles. 
The supercluster area is still a rather polymorphic assembly 
of matter, connected by means of filamentary extensions. 
The most prominent concentrations, the clusters, roughly 
follow these weblike structures. 

Interestingly, we find that the mean mass of all the clus- 
ters in our simulation volume is Mc; = 9.4 x 10^^h~^MQ 
while those residing within the realm of superclusters have 
an average mass Md — 3.6 x lO^*/i~^M0. The fact that su- 
perclusters contain more massive clusters is par tly a result of 
the stronger c lustering of higher mass clusters (|Kaiseilll984l : 
lBahcalj[l988l ). An additional factor is the active dynami- 
cal environment inside superclusters. Because of the high 
concentration of subclumps, these are continuously merging 
and falling into ever more massive clumps. Clusters will be 
centres of action and thus grow rapidly in mass. Meanwhile, 
lower mass and lower density clumps are more liable to lose 
mass or even to get gradually dismantled by the prevailing 
strong tidal forces in and near the superclusters. 



9 SUPERCLUSTER SUBSTRUCTURE AND 
MULTIPLICITY FUNCTION 

The third aspect of the internal evolution of superclusters 
is that of their substructure. Because superclusters are usu- 
ally identified via their cluster content we focus on their 
multiplicity Nsc, ie. the number of clusters they contain. 
Clusters are taken to be virialized subclumps with a mass 
higher than 3 x lO^^/i'^M©, the low-mass threshold for the 
virialized groups at a = 1. We should take into account that 
the definition of the multiplicity Nsc depends on the cluster 
mass threshold in our sample. 

Fig. [15] gives an impression of the substructure of one 
of the superclusters in our sample, at a = 1. The cluster 



9.1 Multiplicity Evolution 

One of the principal findings of our study is that, without 
exception, at a = 100 all superclusters in our sample have at- 
tained a multiplicity one. By that time, they all have evolved 
into compact and smooth density concentrations, akin to the 
one seen in Fig. [3] The hierarchical development of the su- 
percluster involves the gradual merging of its constituent 
subclumps into one condensed object. As a result, we see 
that superclusters which at the present epoch contain sev- 
eral to dozens of clusters ultimately end up as an object 
of unit multiplicity. Fig. [16] reveals the systematic evolution 
towards this configuration. 

Turning to the present epoch, we assess the multiplicity 
function for the 17 most massive superclusters, those whose 
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Figure 16. Supercluster multiplicity. Left panel: the cumulative 
multiplicity distribution isc(Nsc) of superclusters containing 
Ngc or more clusters at four different expansion factors: a = 1, 
a = 2, a = 5 and a = 10. At each timestep the distribution 
concerns the most massive superclusters in the simulation, with 
M> 10^^/i~^Mq. Righthand: Supercluster multiplicity Ngc as a 
function of the mass M of the supercluster. 



mass is in excess of M= 5 x lO^^/i~^M0. When evaluating 
the (cumulative) multiplicity distribution, the number den- 
sity of superclusters with more than Nsc clusters (Fig. 1161 
top-left row), we see that half of the superclusters have 10 
or more members. Another important trend, not entirely 
unexpected, is that larger and more massive superclusters 
contain a higher number of cluster members. One can im- 
mediately infer this when looking at the multiplicity Nsc 
against supercluster mass Msc, plotted in the righthand row 
of Fig. 1161 This trend is particularly strong for superclusters 




Figure 17. Mean mass of the cluster members of superclusters 
in our sample, including (diamonds) and excluding (triangles) the 
central mass concentration. 



with Nsc < 10. For a reason which we do not entirely un- 
derstand, the multiplicity seems to level off for more massive 
and larger supercluster complexes. 

At first gradually, at a = 2, and later more radically we 
see a strong change in the multiplicity Nsc of superclusters 
as more and more clusters within their realm merge and mix 
with the central matter concentration. The lefthand column 
of panels shows the systematic decline of high multiplicity 
objects. At a = 2 we can no longer find superclusters with 
more than 10 cluster members. At that epoch 2 of the 17 su- 
perclusters have collapsed into a single object of multiplicity 
one. Towards the later epochs the multiplicity of superclus- 
ters quickly declines further, and at a = 10 there are no 
superclusters around with more than 5 clusters. 

When looking at the corresponding frames for the mul- 
tiplicity as a function of supercluster mass, we see that the 
number of superclusters of unit multiplicity is continuously 
increasing. Meanwhile, the superclusters that still have more 
than one cluster member appear to have less and less clus- 
ters within their realm. 

An additional interesting issue is that of the masses of 
the clusters that populate superclusters. In Fig. [17] we can 
see that the mean mass of clusters in superclusters is con- 
tinuously growing, from Mcj = 7 x lO^''/i~^M0 at a = 2 to 
Mc; = 1.6 X lQ^^h~^MQ at a = 10. However, this reflects 
the continuous growth of the central mass concentration, ac- 
companied by the decreasing supercluster multiplicity: the 
mean mass is becoming more and more dominated by that 
of the central supercluster core. Meanwhile, the mass of the 
remaining clusters outside the core is decreasing. In other 
words, the first clusters to merge with the supercluster core 
tend to be its most massive companions. The ones that ac- 
crete at a later epoch have a considerably lower mass. This 
is in accordance with our observation in sect. |3] (see Fig. O. 



9.2 Multiplicity Criterion 

Because we have used a physical criterion for the identifi- 
cation of clusters and superclusters, we assume the inferred 
multiplicities to be close to the one found in the observed re- 
ality. To some extent supercluster multiplicity estimates will 
depend on the supercluster identification procedure. There 
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is certainly a dependence on percolation rad ius when identi- 
fied on the basis of a percolation c riterion jZeldovich et alj 
ll982l : IShandarin fc Zeldovichlll989^ ■ 

This may be one of the reasons why our measured su- 
percluster multiplicities differ from those obtained in other 
studies (in addition to the dependence on cluster mass 
threshold). For example. IWrav et al.l l|2006l ) found superclus- 
ters with more than 30 members. This certainly relates to 
the choice of linking length for defining the supercluster: de- 
pendent on the linking length they found maximum super- 
cluster sizes ranging from ~ 30/i~^Mpc to ~ 150/i~^Mpc. 
The latter are much larger than the superclusters we find 
according to our definition based on binding-density. 



10 SHAPLEY-LIKE SUPERCLUSTERS 

The Shapley concentration, first noted by IShaplevI l|l930l ). 
is one of the most outstanding s uperclus t er complexes 
out to z = 0.12 (see Raychaudh ur^ Il989l : lEttoril [l997l : 



iQuintana et al ] |2000l : Froust et al .1 12009 ). It amasses at least 
30 rich Abell galaxy clusters in a core region of ~ 25/i~^Mpc 
and is located at a distance of ~ 140/i~^Mpc. Its total mass 
is estimated to be ~ 5 x l0^^h~^Mp, wi t hin a radius of 
~30/i"^Mpc (see e.g. IProust et all l2006l : iMufioz fc Loebl 
|2008| ). Almost as massiye is another similar assembly of 
massive clusters in the local Universe, the Horol ogium- 
Reticulum supercluster (see e.g. iFleenor et al.|[2005^ . Both 
may have a major influence on the motion of the Local 
Group with respect to the background Universe (s ee e.g. 
IPlionis fc Valdarninlll99ll : lKocevski fc Ebeling||2006l ). 

We investigated in how far we can find back the equiv- 
alents of the Shapley and Horologium superclusters in our 
simulation. On the basis of the mass functions determined 
in sect. 16.21 we estimate that the most massive supercluster 
present in a volume akin to the local Universe {z < 0.1) 
may have a mass of ~ 8 x lO^^/i~^M0. This mass is slightly 
larger than the one in the b ound region of the Sh apley con- 
centration as determined bv lDunner et al.1 l|2008 l). Extrapo- 
lating this conclusion, we find that one may typically find 
two Shapley- like superclusters out to z ~ 0.1, a volume of 
the size of the Local Universe. 

When turning towards the multiplicity of the detected 
simulation superclusters, we find from Fig. [16] that the 
largest supercluster in the Local Universe would have 15 
members. A Shapley-like supercluster would have a radius 
of ~ 14/i~^Mpc and host between 10 to 15 members, close 
to the numb er found in the bou nd region of the Shapley 
supercluster (|Diinner et al.l |2008| ) , which contains ~ 1/3 
of the clusters tra ditionally assigned to this structure (e.g. 
iProust et al.ll200^ ). In the observational reality of our Local 
Universe {z < 0.1), we find 5 superclusters with 10 or more 
members, the large st one containing 12 members (see e.g. 
lEinasto et al]ll994l ). 

This brings us to the issue of the extent to which the 
supercluster population in our LCDM simulation at a = 1 
resembles the one seen in our nearby Universe. To this end 
we are in the process of translating our theoretical superclus- 
ter criterion into one that would be able to identify struc- 
tures along the same lines in large magnitude-limited galaxy 
reds hift surveys. By a pplying this to the 6dF galaxy sam- 
ple n ones et akl |2009| ) and the SDSS DR7 galaxy sample 



(|Abazaiian et al]|2008l ). we plan a detailed comparison of 
morphology, size and spatial distribution of the identified 
superclusters with those seen in our simulation. First indica- 
tions from similar comparison betwee n superclusters in the 
2dFGRS sa mple JCoUess et al.1 [20031 ') and the Millennium 
simulation l|Springel et al.l 2005h . based on a more conven- 
tional supercluster identification process, see ms to indicate 
that t here are no discernable discrepancies (jEinasto et al.l 
l2007bl ). 



11 CONCLUSIONS 

In this work, we have followed the evolution of bound ob- 
jects from the present epoch up to a time in the far future of 
the Universe, at a = 100, in a standard ACDM {ilm,o = 0.3, 
f^A.o = 0.7 and h = 0.7) Universe. We contrasted the ex- 
ternal global evolution of the population of bound objects 
with their vigorous internal evolution, starting from the con- 
tention that in a dark energy dominated Universe they have 
the character of island universes. Within such a Universe 
we expect them to become increasingly isolated objects in 
which cosmic evolution proceeds to the ultimate equilibrium 
configuration of a smooth, spherical, virialized and highly 
concentrated mass clumps. We identify the most massive of 
these objects with superclusters. 

For the external evolution, we investigate the spatial 
distribution and clustering of bound objects and superclus- 
ters, along with the weak change of their mass functions. 
To assess the internal structure, we have looked into their 
rapidly changing shape, their evolving density profile and 
mass concentration and the level of substructure of the su- 
perclusters in terms of their multiplicity, i.e. the number of 
clusters within their realm. 

We defined the bound structures by the density crite- 
rion derived in Paper I (see Eqn. I13|l. and identified them 
from a 500h~^ Mpc cosmological box with 512'^ dark matter 
particles in a ACDM (^2^,0 = 0.3, ^a.o = 0.7 and h = 0.7) 
Universe. We ran the simulation up to a = 100, which is a 
time where structures have stopped forming. We used the 
HOP halo identifier in order to identify independently viri- 
alized structures, at each of the five timesteps we have ana- 
lyzed in detail: a = 1, 2, 5, 10 and 100. 

The main results of the present study can be summa- 
rized as follows: 

• While the large-scale evolution of bound objects and su- 
perclusters comes to a halt as a result of the cosmic accelera- 
tion, their internal evolution continues vigorously until they 
have evolved into single, isolated, almost perfectly spherical, 
highly concentrated, virialized mass clumps. This develop- 
ment is very strong between a = 1 and 10, and continues up 
to a = 100. 

• The marginally bound objects that we study resemble 
the superclusters in the observed Universe. While clusters of 
galaxies are the most massive, fully collapsed and virialized 
objects in the Universe, superclusters are the largest bound 
-but not yet collapsed- structures in the Universe. 

• The superclusters are true island universes: as a result 
of the accelerating expansion of the Universe, no other, more 
massive and larger, structures will be able to form. 

• While the superclusters collapse between a = 1 and 
a — 100, their surroundings change radically. While at the 
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present epoch they are sohdly embedded within the Cos- 
mic Web, by a = 100 they have turned into isolated cosmic 
islands. 

• The large scale distribution of bound objects and super- 
clusters (in comoving space) does not show any significant 
evolution in between a = 1 and a — 100. The cluster and 
supercluster correlation functions do not change over this 
time interval, and retain their near power-law behaviour. 
Superclusters remain significantly stronger clustered than 
the average bound object, with a supercluster correlation 
length of 23 ± 5/i~'^Mpc compared to ro « 11.5/i~^Mpc for 
the fuU bound object distribution in our simulation sample. 

• The mass functions of bound objects and superclusters 
hardly change from a = 1 to a = 100, as we expect on 
theoretical grounds. The mass functions in the simulations 
are generally in good agreement with the theoretical pre- 
dictions of he Press- Schechter, Sheth-Tormen, and Jenkins 
mass functions. At a = 1, the Sheth-Tormen prescription 
provides a better fit. At a = 100, the pure Press-Schechter 
function seems to be marginally better. This may tie in with 
the more anisotropic shape of superclusters at a = 1 in com- 
parison to their peers at a = 100. 

• The change in the internal mass distribution and that 
in the surroundings is directly reflected by the radial density 
profile. Without exception towards a = 100 all objects attain 
a highly concentrated internal matter distribution, with a 
concentration index c = 0.2. In general, the vast majority of 
objects has evolved into a highly concentrated mass clump 
after a = 10. 

• The mass profile in the outer realms of the supercluster 
changes radically from a = 1 to a = 100. At a = 1 it is rather 
irregular, while there are large differences between the indi- 
vidual objects. This is a reflection of the surrounding inho- 
mogeneous mass distribution of the Cosmic Web. In between 
a = 5 and a — 10, nearly all superclusters have developed a 
smooth, regular and steadily declining mass profile. 

• The inner density profile steepens substantially when 
the inner region of the supercluster is still contracting. On 
the other hand, when at a = 1 it has already developed a 
substantial virialized core, the inner density profile hardly 
changes. 

• As a result of their collapse, the shapes of the bound ob- 
jects systematically change from the original triaxial shape 
at a = 1 into an almost perfectly spherical configuration at 
a = 100. For example, at a = 1 their mean axis ratios are 
({s2/si), {s3/si))= (0.69,0.48). At a = 100, they have mean 
axis ratios of ({s2/si), (53/51))= (0.94,0.85). 

• At the current epoch the superclusters still contain a 
substantial amount of substructure. Particularly interesting 
is the amount of cluster-mass virialized objects within its 
realm, expressed in the so called multiplicity function. Re- 
stricting ourselves to superclusters with a mass larger than 
5x 10^^ ft" ^ Mo, of which we have 17 in our simulation sam- 
ple, we find a multiplicity of 5 to 15 at the current epoch. As 
time proceeds there is a systematic evolution towards unit 
multiplicity at a = 100, following the accretion and merging 
of all clusters within the supercluster's realms. 

• In a volume comparable to the Local Universe {z < 0.1) 
we find that the most massive supercluster would have a 
mass of ~ 8x IO^'^/i-^Mq. This is slightly more massive than 
the m ass of the Shapley Supercluster given in lDiinner et al.l 
l|2008l ). When turning towards the multiplicity, we find that 



the largest superclusters in the Local Universe would host 
between 10 and 15 members, close to the nu mber found in 
the bound region of the Shapley supercluster (|Diinner et al.l 
[2008) (which contains ~ 1/3 of the clusters trad itionally 
assigned to this structure, e.g.. IProust et al] (|2006l )). 

While in this study we have addressed a large number 
of issues, our study leaves many related studies for further 
investigation. One of the most pressing issues concerns an 
assessment of the nature of our supercluster objects. This 
involves a comparison with other supercluster definitions, in 
particular in how far our density-based definition fares at 
the earlier epochs when most objects of similar mass will 
have distinct anisotropic shapes. 

Also, while here we follow the evolution of superclusters 
in the standard ACDM Universe, in an accompanying pub- 
lication we will systematically address the influence of dark 
matter and dark energy on the emerging supercluster. There 
we found that dark matter is totally dominant in determin- 
ing the supercluster's evolution. For a preliminary, detailed 
report on the analysis of the role of the cosmological con- 
st ant in the formatio n and evolution of structures, we refer 
to lArava-Meiol (l2008h . 
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excess of the bound sphere and its dimensionless radius f : 



m = ^ m ■ 



(AlO) 



APPENDIX A: LINEAR DENSITY EXCESS OF 
CRITICAL SHELL 

In order to find the hnear density excess So for the critically 
bound shell, we evaluate its evolution at early epochs (a <^ 
1). At these early times - when density perturbations are 
still very small, 5 <^ 1 - the linearly extrapolated density 
excess (Eqn. [TJJ , 



(5(a) = D{a)So 



(Al) 



represents a good approximation for the (real) density of the 
object (Eqn.[l3j. 

At early times the Universe is very close to an Einstein- 
de Sitter Universe, expanding according to 



a{t) 



2/3 



where t, is the characteristic expansion time. 
Following this cosmic evolution. 



to first order yields 
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This approximation, neglecting contributions of order 5^ 
and higher, is reasonably accurate for density perturbations 
5^1. The evolution of the critical shell's radius f(f), the 
solution to the energy equation eq. [S]for iJ = £" = — |, is 
given by the integral expression 

Jo (l-r)VfT2' 



(A5) 



Retaining the two lowest-order terms of this equation. 
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we find that 
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The front factor on the righthand side of this equation 
should be equal to that of Eqn. (|A4[I . so that the charac- 
teristic time tf is found to be 



with that of Eqn. (|A4|I we find 



1 1/2 



Combining expression 
that 



1/2 



3A ^ 20^ 



(A9) 



^2Q,K,o , 

resulting in the following relation between the early density 



APPENDIX B: PRESS-SCHECHTER 
MODELLING OF SUPERCLUSTERS 

According to the P ress-Schechter fo r malism 
(Press fc Schechtd Il974l: iPeacock fc Heavens! I199G| : 
iBond et all Il99ll : IShethI [l998l 'l. the comoving number 
density of halos of mass M at redshift in a cosmic 
background of density pu, is given by the expression 



An 
dM 



pu 



TT M2 cr(M, 



dln(j(M, 2) 



dlnM 



(Bl) 



in which 5c/ quantifies the relative critical overdensity 5c of 
collapse with respect to the variance of density fluctuations 
(j(M, z) on a mass scale M . For a scenario with a power 
spectrum P{k) of the linear density held. 



(B2) 



{M) = 471 / P{k)Lo{kr)k^Ak, 



where uj{kr) is the Fourier-space representation of a real- 
space top-hat filter enclosing a mass M in a radius r at the 
mean density of the Universe, which is given by 



u){kr) 



sin(fcr) cos(fcr) 
(fer)3 (fcr)2 



(B3) 



The critical spherical collapse overdensity value 5c has 
a weak dependenc e on the cosmological background 



a weak dependenc e on the cosmolog i cal background 
(Gnnn fc GottI Il972l: iLacev fc Cold Il993l : lEke et al.l Il996l : 
Kitavam a fc Sutcilll996l ). Useful fitting formulae for the lin- 
ear spherical model col l apse value 5c were obtained by 
iNavarro. Frenk fc Whit3 l|l997h for = FRW Universes 
and for flat Universes, 



0.15(12^)2/3r!^oi85 
0.15(12^)2/3^!^°'' 



<1,Q.A = 0, 



5c{yim) = 

Sim + S7a = 1. 

(B4) 

According to this formula 5c = 1.68 at a = 1 (for Q.m ~ 
0.3), while at a = 100, when Q,,, = 0.43 x 01"^, one finds 
5c — 1.56. A similar expression for the critical virial density 

1998), 



has been given byfervan fc Norman 
1871^ + 82(n„ - 1) - 



- 39(n 

Sim. + S^A 



1871^ + 60(0™ - 1) - 32(0™ - if 



(B5) 

While most applications of the (extended) PS for- 
malism assume perfectly spherical collapse, we know that 
generic gravitational col lapse of primordial den sity peaks 
proceeds anisotropically. ISheth fc TormenI l| 19991 ) improved 
the PS formalism by taking into account the anisotropic col- 
lapse implied by the anisotropic primordial shape of density 
peaks and the anisotropic tidal stresses imparted by exter- 
nal mass concentrations. Mod elling this by means of the el- 
lipsoidal collapse model (e.g.. lLvnden-Belllll964l : Ilckelll973l : 
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White SidI 19791: lEisenstein fc Loebll 19951 : iBond fc MversI 
19961 : iDesiacquea 120081 ) they showed this translates into a 
more fuzzy moving collapse density barrier. The resulting 
mass function, 



1 + 



M2 (7(M) ■ 



dhia(M) 



dlnM 



(B6) 



g2cr^(M) 



with a = 0.707, p — 0.3 and A ~ 0.322, gives a substantially 
better fit to the mass functions obtained in N-body simula- 
tions. In comparison with the standard PS mass function, ST 
predicts a higher abundance of mass ive objects and a sma ller 
number of less massive ones. Later. [Jenkins et al.l (|200lh re- 
ported a small disagreement with respect to N-body sim- 
ulations: underpredictions for the massive halos and over- 
predictions for the less massive halos. They suggested the 
alternative expression: 



(B7) 



dnj ^ pu d In o-(M) (_| ih^-i+bd^ 

dM AP dlnM ^ 
with A = 0.315, B = 0.61 and e = 3.8. Note, however, 
that their expression does not depend explicitly on 5c. They 
showed that "for a range of CDM cosmologies and for a suit- 
able halo definition, the simulated mass function is almost 
independent of epoch, of cosmological parameters, and of 
initial power spectrum". 
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